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Abstract 

We investigate the propagation of random fluctuations through biochemical networks in which 
the concentrations of species are large enough so that the unperturbed problem is well-described by 
ordinary differential equation. We characterize the behavior of variance as fluctuations propagate 
down chains, study the effect of side chains and feedback loops, and investigate the asymptotic 
behavior as one rate constant gets large. We also describe how the ideas can be applied to the study 
of methionine metabolism. 
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1 Introduction. 



There are two different natural contexts in which stochastic dynamics arises in the study of bio- 
chemical reaction networks. In the first, the stochastic chemical dynamics arises from the ran- 
domness inherent in the formation and breaking of chemical bonds. This "intrinsic stochasticity" 
is particularly relevant when the numbers of molecules are small such as in gene transcription 
and small gene regulatory networks where the mean concentrations no longer faithfully model the 
chemical dynamics. There is a large literature in this field beginning with |4|, including [ 12J . dSJ, 
and recently exemplified by |7|| 13|. In this setting, one typically assumes that the reaction system 
is described by a Poisson process that models individual discrete chemical reactions. One then de- 
rives a partial differential equation for the time evolution of concentration densities. As all species 
have their own intrinsic stochasticity, this partial differential equation is parabolic with a uniformly 
elliptic generator. 

In the second context, which is our focus here, one wants to investigate the response of a large 
biochemical system to external excitation. It is natural and theoretically useful to consider stochas- 
tic excitations and to study the emergent properties of the network as the random fluctuations 
propagate through the system. Here the randomness is a tool used to study the out-of-equilibrium 
dynamics of the biochemical system. In this setting, we assume that the concentrations are large 
enough so that the unperturbed dynamics is faithfully modeled by ordinary differential equations. 
Typically, one in interested in perturbing a single (or small number of) input(s) with white noise. 
Hence, the perturbed problem becomes a stochastic differential equation with a hypoelliptic gen- 
erator. 

The central biological goal driving our work is to understand the behavior of biochemical 
systems in cells, which in vivo are exceptionally large and complicated. A metabolite can be the 
substrate for many different enzymes and participate in apparently unrelated reactions. Individual 
reactions usually have nonlinear kinetics catalyzed by enzymes that are themselves inhibited or 
excited by products or distant substrates in the network. Cells and tissues differ because the genes 
that code for certain enzymes have tissue specific expression patterns and biochemical substrates 
themselves also influence gene expression. Further, each cell's environment, its inputs and outputs, 
and its internal state (e.g. stage of cell cycle) are not constant but vary in time. This continual 
variation affects both the concentrations of substrates and the expression of genes that catalyze 
particular reactions. Thus, the gene-biochemical network should not be viewed as a fixed object 
but as one that is continuously changing. 

For each signal, either external or internal, that causes a particular cell to dramatically change 
its operation, there are two natural questions. First, how does the gene-biochemical network re- 
spond to accomplish the change? Second, how does the network enable the cell to maintain home- 
ostasis in all its other operations despite the change? One would like to understand the structural 
and kinetic principles that allow the network to accomplish both tasks simultaneously. We take 
two distinct approaches to this biological goal. First, we study how fluctuations propagate through 
relatively simple systems. We are interested in discovering how different network geometries mag- 
nify or suppress fluctuations since this may give clues to why biochemical networks look the way 
they do. Secondly, we apply fluctuations to in silico representations of specific biological net- 
works. By observing how fluctuations propagate we can identify reactions or subsystems that are 
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buffered against such fluctuations, i.e. are homeostatic. Then, through in silico experimentation 
(e.g. removing particular reactions), we can take the system apart piece by piece to discover the 
regulatory mechanisms that give rise to the homeostasis. 

In this paper, we develop fluctuation theory for chemical reaction systems for which each com- 
plex (in Feinberg's terminology, [5|) consists of a single chemical species and the kinetics are mass 
action. Thus the corresponding differential equations are linear, so we refer to such networks as 
linear SSC (single species complexes) networks. Because of the linearity, the technical difficulties 
involved in studying the associated stochastic processes are minimized. Thus, linear SSC systems 
are an excellent arena for investigating the effects of network geometry on the propagation, mag- 
nification, and suppression of fluctuations. The principles discovered then become the natural goal 
for generalization to nonlinear settings 

To see the kinds of questions we want to ask, consider a simple chain with a side branch. The 
chemical species are Xi, . . . , X„ , Xg ; the corresponding concentrations are denoted by xi , . . . , x„ , x, 
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The chain has a constant input /, which is perturbed by some random process, in this case, white 
noise. If the input is fluctuating, then each of the concentrations will fluctuate as will the fluxes, 
kiXi. Suppose the side chain is absent. Then, will the variations of the fluxes increase, decrease, 
or stay the same as we move down the chain? Does the answer depend on the rate constants kil If 
the side chain is present, does it affect the variances of the fluxes on the chain? If so, what is the 
effect of the size of L. 

The chemical reaction diagram corresponds to set of differential equations for the concentra- 
tions and, similarly, the diagram with stochastic forcing corresponds to a system of stochastic 
differential equations (SDEs): 

dxi = (/ — kiXi)dt + adB{t) 

±2 = kiXi - LX2 - k2X2 + ks,2Xs 

is = k2X2 - ksxs 



These SDEs in turn give rise to a stochastic process on the state space ]R"+ . We prove that this 
stochastic process has a unique stationary measure. Intuitively, this means that at large times the 
joint distribution of values of the concentrations becomes independent of the initial condition and 
independent of time. That is, the statistics converge to an equilibrium distribution. The variances of 
the concentrations referred to above are the variances of the marginal distributions of this measure. 
We prove the existence of the stationary measure for linear SSC systems in Section 2.2. In Section 3 
we study the propagation of fluctuations in chains. In Section 4 we study the effects of side reaction 
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systems, and feedback loops. In Section 5 we ask what happens to variances in the asymptotic limit 
as one of the rate constants goes to oo, corresponding to a very fast reaction. In Section 6 we show 
how to use the fluctuation theory ideas to investigate methionine metabolism. 

It is important to note that our, goals, methods and results are different from those in classical 
biochemical control theory lfTT1l .l3IJ9l.l21 1. In that theory one takes a system at a fixed steady 
state, makes a small perturbation in a parameter (perhaps an input), and allows the system to relax 
to a new steady state. By comparing the new value of a variable (a concentration or flux) to the 
old value, one computes the percentage change of the variable per unit percentage change in the 
parameter. Technically, one is computing a partial derivative. This kind of sensitivity analysis 
gives good information about local, linearized behavior near the initial steady state. By contrast, 
we are concerned with responses to large scale fluctuations in inputs. Technically, this means com- 
puting properties of the distribution of each concentration or flux from properties of the stationary 
measure. 

It is true that the classical biochemical control theory can be made "stochastic" in the following 
way. Suppose that the system has input / and is at steady-state. Consider the same system with 
input 1 + 1], where is a random variable drawn from some density. For each r] we let the system 
relax to steady state and measure the value, v, of some concentration or flux, f is a random variable 
and comparing it's variance to the variance of i] gives information about how much the steady state 
value of V changes as t] changes. However, this modified biochemical control theory often gives 
completely different answers from the fluctuation theory that we are developing and the differences 
are biologically significant. Consider the chain (without the side chain) in the example above. If the 
input is / + ?7, then, at steady state, the flux knXn must equal / + r], so Var{knXn) = Var{ri); thus 
the variance remains constant down the chain. By contrast, we will see below that in our fluctuation 
theory, under a variety of reasonable assumptions, that the variances of the fluxes decrease as one 
proceeds down the chain. This result is interesting from a biological point of view because it says 
that one way to stabilize the flux out of a chain (i.e. small variance) is to have many intervening 
biochemical steps between the input and the output. 

2 SSC networks and the stationary measure 

In this section we introduce the class of chemical reaction systems that we will study and prove 
the existence of a stationary measure. 

2.1 SSC systems with mass actions kinetics 

Throughout we use the terminology introduced by Horn, Jackson, and Feinberg fT(51f51. Let m be 
the number of chemical species. We shall study chemical reaction systems such that each complex 
contains a single chemical species and refer to such systems as SSC networks. In the sequel, we 
use only the statements in Lemma 2.3. 

Lemma 2.1 (Deficiency of SSC networks). An SSC network has deficiency zero. 
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Proof. Suppose the network has a single linkage class and let S denote the stoichiometric subspace. 
Choose any reaction in the network, Xj — > Xj. Here we have two complexes and one reaction 
vector in S. Thus, if there are no other complexes, we are done. Because the diagram has one 
linkage class, if there are other complexes, then there must be one, call it X^, with an arrow to 
or from either Xi or Xj. This adds one complex and one dimension to S since X^ is not a linear 
combination of X^ and Xj. Continuing in this manner until we have exhausted all the complexes, 
we see that the number of complexes is one greater than dim{S}. Since there is one linkage class 
the deficiency of the network is zero. The case where there is more than one linkage class follows 
easily because the reaction vectors in different linkage classes are orthogonal. □ 

We will concentrate on SSC networks containing the zero complex that have one linkage class. 

Lemma 2.2 (Dimension of S). In an SSC system containing the zero complex with one linkage 
class, dim{S} = m. 

Proof. Since the network contains the zero complex, the number of complexes, n, is one greater 
than the number of species, m. If s = dim{S}, then, by Lemma ITTl = 72 — s — l = m — s, so 
s = m. □ 

We assume mass action kinetics so the differential equations governing the system are linear: 

x{t) = Ax{t) + J, (1) 

where A G M™^™ and x{t), I E M™. The matrix A is the matrix of rate constants for the system 
and the vector / represents any constant flow into the species of the system from the zero complex. 
Thus the components of / are non-negative. We denote the open positive orthant and its closure by 
W^Q and M>o, respectively. 

Lemma 2.3. If a linear SSC system is weakly reversible and contains the zero complex, then 

(a) The differential equations O have a unique equilibrium which is globally asymptotically 
stable and contained in IR>o- 

(b) The eigenvalues of the matrix of rate constants. A, have strictly negative real parts. 

(c) For all vectors v G ]R>o, we have e^*w ■ Cj > 0. 



Proof. Part (a) is a special case of the zero deficiency theorem Q. Since A is the Jacobian at the 
equilibrium point, (b) follows from (a) and linearity, (c) holds because R>o is invariant under the 
flow of the differential equation. □ 
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2.2 The Stationary Measure. 

Consider the following weakly reversible SSC system with mass action kinetics, input vector /, and 
matrix of rate constants A perturbed by a mean zero, finite variance stationary stochastic process 

x{t)=Axit) + I + m , ^2) 
x(0) = Xq ■ 

From this definition and the stationarity of ^ (t) one easily sees that Q generates a time-homogeneous 
Markov process. 

Theorem 2.4. The process x*{t) = x*{t, ^) defined by 

x*{t,0= [ e^^'-'^Ids+ f e^'-'-'^^.ds (3) 

J —oo J —oo 

is a stationary solution to ©. Furthermore given any initial condition Xq, ifx{t, Xq,^) is a solution 
to equation © then x(t, Xq, converges to x*{t,^) as t oo in that 

E\x{t,xo,0 -x*it,0\^ ^0 as t-^oo. 



Proof. Observe that for any t, r G M, 

J —oo J ~oo 

J —oo J —oo 

This can be written succinctly as 

{erX*){t,i)=x\t,erO (4) 

where the shift 6t is defined by {6tf){s) = f{t + s) for all s, t G M and functions / on R. Hence 
for any ti < ■ ■ ■ <tn, 

(x*(r + ti,0,--- ,a;*(r + t„,0) = (a;*(ti, ^,0, " ■ ■ , x*(t„, ^.0) • 

Since ^ is a stationary process, the distribution of the right hand side is independent of r which 
proves that x* is stationary. Clearly, x*{t,^) is a solution in that x{t, x*(0, ^) = ^). 

We now turn to convergence. It follows from Lemma 2.3(b) that there are constants a, M > 
such that ||e^*|| < Me""* for all t > 0. Subtracting the solution of (2), 

x{t,xo,0 = e'^'xo+ [ e^^'-'Hds+ [ e^^'-'^^sds, (5) 
Jo Jo 



from x*{t), squaring, and taking expected values gives, 



E\x{t,xo,0 - x*it^O\ <3||e^*||^|xor + 3E 



<3M"|xo|V2"* + 



+ 3E 



<3M2|a;o|V2"* + 



3M^ J 



oo 

21 T\2 



+ 3E 



J — oo 



oo 
21 r|2 



3M^|/| 



Thus, E|x(t,xo,0 ^ as t ^ oo . 



□ 



Remark. If one takes expectations on both sides of equation (5), one sees immediately that the 
model is consistent in the mean, that is, the mean of the perturbed problem is equal to the solution 
of the unperturbed problem. 

If instead of random perturbations given by the vector we had allowed the system to be 
perturbed by independent white noise processes, we arrive at the following system of Ito stochastic 
differential equations: 

[ dx{t) = {Ax{t) +I)dt + T.dB{t) , 

where S G M™^^ and B{t) is standard p-dimensional Brownian motion. The following theorem is 
proved in the same manner as Theorem 12. 41 



Theorem 2.5. The process x*{t) = x*(t, B) defined by 

x*{t,B) = f e^^'~'Hds+ f 



T.dB{s) 



(7) 



is a stationary solution to ©. Furthermore given any xq, if x{t, xq, B) is a solution to equation 
^ then x(t, xo, B) converges to x*{t, B) as t ^ oo in that 



E\x{t,Xo,B) -x*{t,B) 



as t ^ oo . 



Proof. The proof is identical to that of Theorem 12. 41 except that the Ito Isometry is used to control 
the expected value of the square of the Ito integral term. □ 

Since x*{t) is stationary, the distribution of x*{t) is independent of t and invariant under the 
dynamics of Q (or More precisely, defining the measure /i(v4) = P(a;*(0) G A) for all 
measurable A C M'", we see that 




F{x{t,y,OeA)fi{dy). 
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Furthermore, /i characterizes the longtime behavior of the solution in that the distribution of 
x{t,Xo,^) converges to as t oo. This follows from E|x(t, Xo,0 ~ and the 

fact that ij{A) = ¥{x*{t) e A) for all t. 

Thus fj. contains information about the average, long-term behavior of fluxes and concentra- 
tions. It will be fi, therefore, which we shall probe in order to gain an understanding of how differ- 
ent graphical structures and asymptotic limits of biochemical reaction systems increase, decrease, 
and otherwise modify the exogenous fluctuations of biochemical reaction systems. Throughout 
the rest of this paper, it is understood that each mean or variance is computed with respect to this 
stationary measure. 

2.3 A General Bound 

We can now prove a simple general bound for the variance of the concentration of any species in 
an SSC system in terms of the variance of the input fluctuations. We assume that the fluctuations, 
^t, are one-dimensional, stationary, mean zero, and finite variance. By taking the expected value 
in equation Q and using that has mean zero one sees that 



is the mean of the i species. 

Theorem 2.6. Let x*(t) be the stationary solution of an SSC system with one input, I, to a single 
species, Xi, that is perturbed by a stationary stochastic process, C,t, with finite variance and mean 
zero. Then for each i. 




(8) 




Proof. Using Lemma 2.3(c) and the Cauchy-Schwarz inequality gives 




The strictness of the inequality follows because is not a constant. 



□ 
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This simple result is all that we need in this paper. An analogous proof works in the more 
general case where there are inputs to more than one species and any number of the inputs undergo 
independent fluctuations. 



3 Chains 

In this section we consider non-reversible chains with mass action kinetics: 

> Xi > X2 >■ • • • > Xm-1 > Xm ^ 0. 

Theorem 12 . 61 allows us to see that variances of the fluxes of the stationary solution decrease as one 
proceeds down the chain. 

Theorem 3.1. Let the input, I, of a non-reversible chain with mass action kinetics be perturbed 
by a stationary stochastic process, ^t, with finite variance and mean zero. Let x* (t) denote the 
stationary solution for the chain. Then, for all i, Var{kiX*) < Var{C,) and 

Var{ki+ix*^-^) < Var{kiX*). (10) 



Proof. From the remark following Theorem 2.4, we know that the mean, m,, of x*(t) is the equi- 
librium value of Xi for the unperturbed problem. For the chain this implies that rrii = p, so the 
bound Var{kiX*) < Var{^) follows immediately from Theorem l2.6l To prove (fTOI) note that the 
input to X2 is 

kixl{t) = I + (A;ix*(t) -/) 

and kixl{t) — / is a stationary stochastic process of mean zero and finite variance. Thus, by 
Theorem 1221 

Var{k2xl) < Var{kix\ — I) = Var{kix\) . 
The input to X3 is k2xl{t), so repeating this argument down the chain proves (fTOb . □ 

Note that the variances of the fluxes are strictly decreasing as one moves down the chain even 
though the means of the fluxes remain unchanged (i.e., equal to /). The next natural question 
is how much do the variances decrease down the chain? This cannot be answered without more 
detailed information about To investigate it, we will perturb the input / by white noise, adB{t), 
which will allow us to use the Ito calculus. 

Theorem 3.2. Let x*{t) be the stationary solution of the linear chain © where the input is per- 
turbed by white noise. We assume that the rate constants, ki, are distinct. Then 

i i ^ 

Var{x*) = ^^PijPir , , , , (H) 
j=l r=l '^^ + '^^ 
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where 



Pij 







^ > 3 
i < j 



(12) 



Proof. The matrix of rate constants, A, is given by 

' -ki 



A 










h 1 —k 



Let P = {pij}. A straightforward calculation shows that the jth column of P is the eigenvector 
of A corresponding to eigenvalue —kj. Thus, D = P^^AP is diagonal. In addition, P takes the 
vector (1, 1, ■ ■ ■ , 1)^ to the vector (1, 0, ■ ■ ■ , 0)^. Using these facts, the formula Q) for x*{t), and 
the Ito Isometry, 

2 



Var{x*) = cT^E e^^'-'^d ■ ddB^j 



a 







■ Cj) ds 


( 


' 1 " 


\ 


2 


P(,D{t-s) 






ds 


\ 


1 


) 





a 



a 



a 



ds 



i i ^ 
X] PijPir u,l. ■ 

Ay-i rvf 



ds 



j=l r=l 



□ 



We assumed that the fc/s were distinct so that the explicit formulas above make sense. It can 
be shown that the variances of the concentrations are continuous functions of the rate constants. 
This fact, together with the bound given by (flOb allows us to conclude that formula (fTTT) has finite 
limits as various subsets of the /cj's become identical. 

We can use the explicit formula dTTl) to answer several natural questions: 
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Example 3.3. (Magnitude of decrease) Theorem 13.11 shows that variances of fluxes are strictly 
decreasing as one moves down a chain. To investigate how much they decrease, consider the chain 
® where m = 2 and the input is perturbed by white noise. Using dTTT) we see that Var(kixl) = 

^ and Var{k2xl) = Thus, 

Var{k2xl) k2 
Var{kixl) ki + k2 

This simple example shows that the ratio of successive variances can be any number between zero 
and one. 

Example 3.4. (Long chains) Assume that ki = k for some fixed A; > and all i. Taking the limit 
of dTTT) is difficult. Instead, since all the ki are equal, an induction proof shows that 

xm = ^- + aj^_ J\t - sy-'e-'^'-^UB{s). 

Using the Ito Isometry, it follows that 

Var (X*) = ,^^(^1^1 

and using Stirling's formula 

Varikx*) ~ a'-^^ + 0{z--'/'). 

Thus the variances decrease to zero in a regular fashion if all of the rate constants are the same. 

Example 3.5. (A small rate constant) Suppose that one rate constant, ki, in a chain is very small. 
Using the explicit formula dTTT) . one can easily compute that 

Var{kiXi) a^^ki + 0{k^), as ki 0, 

Var{kjX*) ~ cr^ifcj + O^kf), as ki 0, for j > i. 
2 

Notice that the small rate constant has the effect of significantly decreasing the variances of the ith 
and all subsequent fluxes while the means of the fluxes remain unchanged. Therefore a small rate 
constant is not "rate limiting" but instead is "variance limiting." 

Example 3.6. (A large rate constant) Suppose that one rate constant, ki, in a chain is very large. 
Again, using (fTTT) . one can compute that 

Var{kiX*) Var{ki^iXi_i), as ki — > oo. 

Furthermore, for all j > i, 

Var{kjX*) Var{kjX*), as ki oo, 
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/ + adB{t) ki k2 h-i h+i 

^ Xi ^ X2 ^ ■ ■ ■ Xi^i ^ Xi^i ^ 

where Xj is from the process arising from the following system: 

This shows that in the asymptotic limit where ^ 00 one can replace the original chain by 
the chain with the substrate X,; removed. Here we implicitly use the fact that since the kinetics are 
linear and hence the concentrations are Gaussian the statistics are determined by the means and 
variances. 



4 Side Reaction Systems and Feedback Loops. 

A side reaction system on a chain is any SSC system that gets its input from a species on the chain 
and has output that flows back into the same species; see Fi gure l4~T] below. 



m ^ k. 



Xo 



k. 



A;, 



Side Reaction System 



Figure 4.1: A side reaction on a linear chain 



Note that there must be a species within the side reaction system whose output flows to Xi with 
some rate constant, ^3. Define Y to be that species. The SDE governing the behavior of Xi{t) is 
then given by 

^xi(t) = / - fcixi(t) - k^xiit) + k^yit) + i(t). (13) 

If Xi is the solution to the above system when there is no side reaction system (i.e. k2 = k^ = 0), 
then 

^Mt) = I-k^Xl{t)+m■ (14) 
at 

Theorem 4.1 (Side reactions lower variance). Let x* and xl be the first components of the sta- 
tionary solutions to (O and (O, respectively, where C,(t) is a finite variance, mean zero, random 
process or white noise. Then, 

Var{kix\) < Var{kix\). 



Proof. We give the proof in the case where ^{t) = adB{t) is white noise; the proof in the general 
case is similar but more complicated [1|. Let z{t) = E,(kiXi{t) — J)^ and z{t) = E,{kiXi{t) — J)^, 
where xi{t) and xi{t) are solutions of (fT3l) and (fT4l) . By theorem |231 and z{t) converge to 
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Var{kixl) and Var{kixl), respectively. We will prove the theorem by comparing the differential 
equations satisfied by z{t) and z(t). 

By using Kolmogorov's backward equation ifTHll . we see that z{t) satisfies 

z'{t) = -2kiz{t) + kfa^. (15) 

Therefore, Var{kixl) = kia'^ /2 because Var{kixl) is the equilibrium value of ( [T5t . Similarly, 
z(t) satisfies 

z'{t) = -2kiz{t) + kla'^ + 2A;iE(A;iXi(t) - I)ik3y{t) - k2Xi{t)), 
and so, by Theorem |231 

Var{kixl) = ^ + ^Hhxl - I){ksy* - k2xl). 

Thus, to complete the proof we need only show that 'K{kixl — I){k^y* — k2xl) < 0. The remark 
following Theorem 2.4 implies that Exl = I/ki and Ek^y* = EA;2X*. Therefore, E/cay* = 
and 

E(fcixl - I){ksy* - k2xl) = ^E {k2xlksy* - kjxf) . 

K2 

By Theorem IT6l E (k^y*)^ < E {k2xlf, so 



\E{k2xlk3y*)\ < {Eklxfy (E kly*^y 
< E klxf. 

Thus, E{kixl - I){ksy* - k2x\) < 0, as desired. □ 



A feedback loop on a chain is an SSC system together with an input from one species on the 
chain, X„, and an output to an earlier species, Xi, see Figure l4!2l 




Subsystem 



Figure 4.2: A chain with a feedback loop 

Theorem 4.2. Let x{t) be the vector of species concentrations for the chain ^ and let x{t) be 
the vector of species concentrations for the chain with feedback loop (Fisure W^ . where ^{t) is a 
finite variance, mean zero, random process or white noise. Then, 

Var{knX*J < Var{knX*J. 
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Proof. Let {V^} be the substrates and B be the matrix of rate constants of the SSC subsystem in 
Figure We suppose that Vj is the species which gives input to Xi with rate constant a. Then 
the input to Xi from the feedback loop is 

f\{t) = ae^^iQ) ■ ej + ac [ e-^(*"'^x„(s) ■ ej ds, 



which depends explicitly only on Xn. If we let R{t) = A;„_ix„_i(t) then the differential equation 

for Xn{t) is Xn(t) = R{t) - CXn{t) - knXn{t). 



Xi — 



kn—l __ i^n 



Xri 



kn 



Subsystem 




n—l 



Figure 4.3: A chain with a side reaction system 



Consider the chain with side reaction system given in Figure 14.31 where the subsystem is 
the same as in Figure l4~2l and the flux to Yi comes from Vj with rate constant a. Let Q(t) = 
kn-iXn-i{t) and P(t) = kn-iVn-iit) be the inputs to X„ in Figure l431 Since the input to the Y- 
chain is fi(t) and the rate constants for the two chains are the same, R{t) = Q{t) + P{t) because 
the differential equations are linear. Thus, the differential equation governing Xn{t) in Figure l4~2l 
is the same as the differential equation governing Xn{t) in Figure l431 Since the system in Figure 
I4.3l is a chain with a side reaction system, the result follows from Theorem 14. II □ 



5 One large rate constant in a general SSC system 

We now consider a general weakly reversible SSC system with input and characterize the effect of 
a large rate constant. 

Theorem 5.1. Suppose that independent white noise processes perturb the inputs to a weakly 
reversible SSC system with m substrates. Let Xa be a particular substrate and suppose that the 
rate constant L for one flux out of Xa to another complex (possibly the zero complex) is large. 
Then, 

Var{xl)r^o(\] asL^oo. (16) 
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Proof. We will assume that one of the perturbed inputs goes directly to X^. The proof of the 
general case is similar. The stochastic differential equation governing Xa{t) is given by 

dxa{t) = [c + Y1 ^i^ii'^) - + K)xa{t)^ dt + adB{t) , (17) 

where L + i^T > is equal to the sum of all the rate constants for reactions leaving X^, C > is 
the input flow to X^ from the zero complex, cr > 0, and q > is the rate constant associated with 
the reaction X, — *• Xq. Solving (flTt for a;* in terms of the x* and using the Ito Isometry, one can 
easily bound Var{xl), 

for some constant p. To complete the proof we will show that Var{x*) < 0{L). 

Let A be the matrix of rate constants for the SSC system. Using the formula dT]) for the station- 
ary solution and the Ito Isometry, one easily calculates: 

\/ar(x*) = a^ ! (e^^^'^^e ■ 6^)^ ds, (18) 



for some vector e. By Lemma lOt b) we know that the real parts of the eigenvalues of A, {Aj}, are 
strictly negative; let A = inf { | Aj | }. There exist positive constants c and M so that for all t — s > 0, 
we have ||e^*^*^*^|| < ce~*^^*^*~''\ Using this inequality in (flSt . we have 



Var{x*) < 



2M A 

In Appendix A we prove that A > 0(1/L), so Var{x*) < 0{L), which concludes the proof. □ 

Example 6.1 (A side chain with a large rate constant) To illustrate the theorem, we consider the 
linear chain with a side reaction given in the diagram in the Introduction. As the rate constant L 
becomes large. Theorem 15. II tells us that Var(x2) < 0(1/ L). Therefore the flux out of X2 down 
the chain has variance Var{k2X*2) < 0(1/ L). By Theorem 12.61 

Var{kiX*) < Var{k2X*2) < 0(1/ L) for aU i > 2. 

Thus, for all i > 2, the means of the fluxes remain equal to /, while the variances of the fluxes go 
to zero as L ^ 00. 



6 Application to Methionine Metabolism. 

The actual biochemical systems involved in cell metabolism are much more complicated and more 
difficult to analyze than the single species systems considered in the previous sections. Consider, 
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for example the diagram in Figure 6.1 that shows the methionine cycle and part of the folate 
cycle. Firstly, most reactions have two or more substrates and many enzymes are inhibited by the 
products of the reactions they catalyze. Thus, the kinetics will be highly nonlinear. Secondly, many 
reactions are catalyzed by two different enzymes that have very different properties. Thirdly, some 
substrates inhibit or activate distant enzymes in the reaction diagram (red arrows in the diagram). 
These long-range interactions make it virtually impossible to intuit the emergent properties of the 
network by tracing influences from point to point. 




cystathionine 



Figure 6.1. Methionine Metabolism. Substrates of the methionine cycle and (part of) the folate cy- 
cle are shown in green and red rectangles, respectively. Enzyme acronyms are in ellipses. Long-range 
interactions are shown by red curves with the arrow indicating activation and the bars indicating inhibi- 
tion. SAM, s-adenosyl-methionine, activates CBS and inhibits BHMT and MTHFR, while 5mTHF, 5- 
methyltetrahydrofolate, inhibits GNMT. 

Epidemiological evidence correlates changes in folate and methionine metabolism to serious 
human health consequences (cancer, heart disease, depression) and there are several important 
public health issues involved in folate supplementation as currently practiced in the United States 
and Canada. Thus, this part of cell metabolism has been the object of numerous experimental 
studies and several modeling studies [fT4lli20ll lfT6lllfT9lllfT7l . Our purpose here is simply to illustrate 
how fluctuation analysis can be used to understand such a complex system. 

The velocities of the individual reactions in the methionine cycle |T7^ are typically highly non- 
linear functions that depend on the concentrations of several substrates. For example, the velocity 
of the GNMT reaction, Vcnmt, depends on SAM, on SAH because of product inhibition, and on 
hmTHF because of a long-range interaction. Because of the complexity and the nonlinearities, 
a rigorous mathematical analysis of this system is beyond current mathematical techniques. Even 
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proving the existence and uniqueness of a stationary measure is a delicate issue. Nevertheless, 
we have investigated this question by numerical computation in the case where the methionine 
input {MET in) is an Ornstein-Uhlenbeck process with mean 100 /xM/hr and standard deviation 
30 (variance = 900). We found that the joint distribution of the substrates does indeed stabilize as 
time gets large, and thus for each concentration or flux, X, we can compute the ratio 

Variance (X) 
^ ~ Variance (AfETm) ' 

which tells us how much X varies compared to the variance of the input. Table 1 shows the 
values of r for two substrates and two fluxes in the case where all the long-range interactions 
are present (regulated) and the case where the long-range interaction are absent (unregulated). 
Methionine is quite stable in both cases but is more stable in the regulated case. SAM is much 
less stable than methionine, which agrees with what is seen experimentally. Notice that in the 
unregulated case, the variances of Vgnmt and Vdnmt are similar, but in the regulated case the 
variance of Vgnmt doubles and the variance of Vdnmt becomes exceptionally small. There are 
good biological reasons why one would want the DNA methylation rate to be stabilized against 
fluctuations in methionine input. Thus, fluctuation analysis shows that this stabilization is achieved 
by the long-range interactions. We have also computed the values or r in all the intermediate cases 
where some but not all of the long-range interaction are present and this has enabled us to quantify 
each of their effects and propose an evolutionary scenario ifTTI . 



Table 1. Values of r 



Methionine 


r 


regulated 


.064 


unregulated 


.082 


SAM 


r 


regulated 


.22 


unregulated 


1.23 


Vgnmt 


r 


regulated 


.15 


unregulated 


0.079 


Vdnmt 


r 


regulated 


0.007 


unregulated 


0.09 



In liver cells the reaction from methionine to SAM is catalyzed by two isoforms of the same 
enzyme, MAT-I and MAT-III, that have very different properties O. MAT-I is inhibited by 
SAM and MAT— III is activated by SAM, and it has been proposed that it is this unusual combi- 
nation that stabilizes the methionine concentration. To test this, we recomputed r after eliminating 
the MAT— III reaction and raising the Vmax of the MAT— I reaction so that it carried the same 
flux previously carried by both. The values in Table 2 show conclusively that, indeed, the presence 
of the MAT-III reaction somewhat destabilizes SAM but greatly increases the stability of the 
methionine concentration. 
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Table 2. Values ofr 



Methionine 


r 


with MAT-III 


0.06 


no MAT-III 


0.16 


SAM 


r 


with MAT-III 


0.22 


no MAT-III 


0.17 



7 Discussion. 

In Sections 2-5, we developed the theory of propagation of fluctuations for the special case of linear 
SSC networks and proved theorems relating variances to network structure. Variances decrease 
down a chain and the presence of side reactions and feedback loops always lowers the variances 
further down the chain. These results are very general in that they hold independent of the choice 
of rate constants. It is tempting to speculate that biochemical systems evolved to be as complicated 
as they are partly because of the homeostasis of exit fluxes achieved by having many intermediate 
steps. We also showed how the large size of a single rate constant affects variances. It is known 
that most of these results generalize to non-linear SSC networks with restrictions on the nature of 
the nonlinearity ||2|. It remains to be seen whether they generalize to networks in which complexes 
contain more than one species. In these highly non-linear contexts, a fundamental mathematical 
issue is the proof of the existence of a stationary measure. 

A reasonable concern with the idealized models in Sections 2-5 is that, under the influence 
of the fluctuations, the concentrations can become negative. By modifying the forcing processes 
appropriately this could have been avoided. However, this would complicate the analysis and 
prevent us from obtaining explicit formulae and straightforward bounds. Since our goal with these 
idealized models is to build intuition and develop general principles, we have purposely avoided 
complicating the analysis. 

In Section 6 we showed how the ideas of fluctuation theory could be used to investigate a 
network of biological interest, the methionine cycle. It is reasonable to ask whether methionine 
input actually fluctuates randomly and if so what are the properties of the fluctuations. There 
are really two answers. The input to the methionine pool in liver cells is certainly continually 
varying. There are large deviations on the time scale of hours depending on the times and content 
of meals. Methionine is always being used for protein synthesis and is being made available by 
protein catabolism, two processes that are themselves variable and not always in balance. The 
methionine available for input to the methionine cycle is also affected by the use of methionine in 
other metabolic reactions. Finally, all these processes are affected by the time- varying regulation 
of the genes that code for the various enzymes. Thus, the first answer is that we don't know 
how methionine input varies but it certainly fluctuates with standard deviations of the order of 30- 
50 /iM/hr on the time scale of hours and with smaller standard deviations on the time scales of 
minutes and seconds. The second answer is that it doesn't matter. We are using the fluctuations in 
methionine input as a probe of the dynamical properties of the system away from equilibrium. Of 
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course, we need to be sure that the properties we find do not depend on the detailed properties of 
the noise. 

For simplicity of exposition, we have discussed the special case where a single input to a bio- 
chemical system is varied. The same ideas can be used to introduce fluctuations in a concentration, 
a flux, or in several places, and then study how the fluctuations propagate throughout the system. 
Understanding the consequences of fluctuations in kinetic parameters is also important because 
kinetic parameters depend on enzyme concentrations and other properties that are variable and 
themselves dependent on time-varying genetic regulation. Analyzing this case requires some tech- 
nical extensions of this work. 

In the Introduction we referred to "intrinsic stochasticity" in contrast to the external stochastic 
forcing that we consider. It would be interesting to consider models with both forms of stochastic- 
ity, and indeed both surely arise in gene networks. In gene networks that are coupled to biochemical 
networks, the intrinsic stochasticity at the gene and gene regulation level can be viewed as external 
stochastic forcing to the biochemical level. Therefore, both types of questions and analyses will be 
necessary to gain full understanding of real biological networks. 

Appendix A 

We derive the bound used in Theorem l5.ll There are two cases which need consideration: 

1. The flux out of Xa with rate constant L goes to another species. This case is handled in 
Theorem A. 1 below. 

2. The flux out of Xa with rate constant L leaves the system. The proof of the result in this case 
is similar to the proof of the theorem below and so the details are omitted. 

Theorem A.l Let A = {aij} be an n x n matrix with the following properties: 

(1) For each i, an < and |ajj| > Yl'jj^i 

(2) ail = —L + au and = L + a2ifor some au < and a2i G R. 

(3) For every L > 0, the real parts of the eigenvalues of A are all strictly negative. 
Denote the eigenvalues of A by {Aj} and let X = inf {|_Re(Ai)|}. Then 

A > 0{1/L), asL^oo. 



Proof. Let B = jA. The eigenvalues of B are {j^Ci : is an eigenvalue of A}. We will use 
the characteristic polynomials of A and B to show that the magnitude of the real parts of the 
eigenvalues of B are no smaller than 0(1/ L^), which implies our result. 

Because L only appears in the first column of A, all 0(1) terms of B occur in the first column. 
Expanding the determinant of B by cofactor expansion along the first column then shows that 
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det{B) must be of order 0(1/L") or 0(1/L"~^). Similarly, the cofactors of B must be of order 
0(1/L"~^) or 0(1/ L"-'"^). Therefore, computing the inverse of B (which exists by assumption 
(3) above) by cofactors, we see that the possible order of the entries of B^^ are 1, L, and L^. 
Therefore, \\B-^\\ < 0{L^). 

One may view 5 as a 1/L matrix perturbation of the matrix C = {q-,}, where cn = — 1, 
C21 = 1, and Cij = for all other entries. Therefore, each eigenvalue, p, of B is an analytic 
functions of 

P = Po + tPi + T^P2 + o(^] , (19) 



where pq is —1 or 0. If po = ~1 there is nothing to prove; so we assume po = 0. If pi = p2 = 
then p = 0{1/L^). However, this would imply that 0(l/p) = 0{L^). Since 1/p is an eigenvalue 
of B^^, this would contradict the norm bound for B^^, above. Thus pi and p2 can not both be zero. 
It remains to be shown that the leading order term in equation ST% can not be purely imaginary. 
We will do this through asymptotic matching. 

Consider two different formulations for the characteristic polynomial of A, pa{x)'- 

p^(x) = det{xln - A) (20) 
= x" + Lm(x) +f(x) (21) 



= x" + ci,„_iLa;" ^ + cq^u-ix^ ^ H h ci^2Lx^ + co,2a;^ 

+ ci^iLx + co,ix + ciflL + co,o, 



(22) 



where u{x) and v{x) are polynomials of degree n — \ that are independent of L, and Cij E M for 
z = 1, 2 and j = 1, ..n — 1 (i gives the power of L and j gives the power of x for the term cijUx^). 
We note that we can not have ci = Co,o = 0, for then there would be a zero eigenvalue, which 
would contradict assumption (3). 

To show that the leading order term in equation (fT9l) is not purely imaginary we will consider 
two cases: pi = and pi 7^ 0. We begin by supposing pi = and p2 7^ 0. Then p = 0{1/ L?') and 
there is a solution to (1^ which is O {1/ L). Putting x = P2/ L into (E3) and setting the equation 
equal to zero gives us: 

^ / 1 V Ci,2pi , Co,2P2 , , Co,lP2 , r , n 

\u) ^^'^^^ " 

Matching like terms in L tells us that ci,o = 0, co,o 7^ 0, and ci,i 7^ 0. Solving for p2 gives us 
P2 = — co,o/ci,i G M. Therefore, p2 has a nonzero real part. 

We now suppose that pi 7^ 0. Because finding an 0(1/L) solution to equation (fT9t is equivalent 
to finding an 0(1) solution to (l2n) . pi must satisfy u(pi) = 0. Let D{x) = xin — A. Then 
u(a;) = D(s)ii + D(x)2i, where D{x)ij is the cofactor of D{x). D{x)ii and D{x)2i differ 
only in the first row, so we may combine the determinants by adding the first two rows. We 
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conclude that 



u{x) — 



—0,22 ~ 0,12 + X — Oi3 — O23 — Oi4 — O24 
-032 -033 + X -034 

—042 —043 —044 + X 



"fln2 



"fln3 



— Oin — a2n 
— Cl4n 



Solving u{x) — for non-zero solutions is therefore equivalent to finding the non-zero eigenvalues 
of the matrix 



A 



022 + O12 O13 + O23 Oi4 -|- 024 
O32 O33 O34 

O42 O43 O44 



fln2 



Cln4 



Oln + 02n 

a4n 



By assumption (1), the diagonal entries of A are non-positive and have magnitudes that are greater 
than or equal to the sums of the magnitudes of all the other terms in that column. Therefore, 
Gershgorin's Theorem says that the non-zero eigenvalues of A, and hence the non-zero solutions 
of u{x) — 0, have strictly negative real part. Thus, Re{pi) ^ 0. This completes the proof. □ 

If the flux out of Xa with rate constant L leaves the system, the only change in the statement of 
the above theorem is that 021 is independent of L. The proof is identical except that u{x) — D{x)\\ 
and so we no longer have to add two determinants together to simplify u{x). 
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